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Abstract 

We introduce a new methodology for a fast and reliable discrimination between ordered and chaotic 
orbits in multidimensional Hamiltonian systems which we call the Linear Dependence Index (LDI). The 
new method is based on the recently introduced theory of the Generalized Alignment Indices (GALI). 
LDI takes advantage of the linear dependence (or independence) of deviation vectors of a given orbit, 
using the method of Singular Value Decomposition at every time step of the integration. We show that 
the LDI produces estimates which numerically coincide with those of the GALI method for the same 
number of m deviation vectors, while its main advantage is that it requires considerable less CPU time 
than GALI especially in Hamiltonian systems of many degrees of freedom. 



1 Introduction 

The fast and reliable discrimination of chaotic and ordered orbits of conservative dynamical systems is of 
crucial interest in many problems of nonlinear science. By the term conservative we characterize here systems 
which preserve phase space volume (or some other positive function of the phase space variables) during time 
evolution. Important examples in this class are N - degree of freedom (dof) Hamiltonian systems and 2N - 
dimensional symplectic maps. As is well known, in such systems chaotic and regular orbits are distributed 
in phase space in very complicated ways, which often makes it very difficult to distinguish between them. 
In recent years, several methods have been developed and applied to various problems of physical interest 
in an effort to distinguish efficiently between ordered and chaotic dynamics. Their discrimination abilities 
and overall performance, however, varies significantly, making some of them more preferable than others in 
certain situations. 

One of the most common approaches is to extract information about the nature of a given orbit from the 
dynamics of small deviations, evaluating the maximal Lyapunov Characteristic Exponent (LCE) o\. If o~i > 
the orbit is characterized as chaotic. The theory of Lyapunov exponents was first applied to characterize 
chaotic orbits by Oseledec jwhile the connection between Lyapunov exponents and exponential divergence 
of nearby orbits was given in S 21 1. Benettin et al. @ studied the problem of the computation of all LCEs 
theoretically and proposed in [7| an algorithm for their efficient numerical computation. In particular, o~\ is 
computed as the limit for t — ► oo of the quantity: 

= 7 ln|^|,i.e.«7x= lirnMt), (1) 

t \\W(0)\\ t-KX) 

where u>(0), w(t) are deviation vectors from a given orbit, at times t = and t > respectively. It has been 
shown that the above limit is finite, independent of the choice of the metric for the phase space and converges 
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to o i for almost all initial vectors w(0) 0; @; 0|- Similarly, all other LCEs, 02, <73_etc. are computed as 
limits for t — > 00 of some appropriate quantities, £2(^)1 -^(i) etc. (see for example [7j). We note here that 
throughout the paper, whenever we need to compute the values of the maximal LCE or of several LCEs we 
apply respectively the algorithms proposed by Benettin et al. 

Over the years, several variants of this approach have been introduced to distinguish between order and 
chaos such as: The Fast Lyapunov Indicator (FLI) 0; [lj; the Mean Exponential Growth of 

Nearby Orbits (MEGNO) [llEi], the Smaller Alignment Index (SALI) [3 H! Hi] , thc Rclative Lyapunov 
Indicator (RLI) (23|, as well as methods based on the study of power spectra of deviation vectors [23 and 
spectra of quantities related to these vectors 

Recently, the SALI method was generalized to yield a much more comprehensive approach to study 



chaos and order in 2N - dimensional conservative systems, called the GALI m indices [27|; |2j. These indices 
represent the volume elements formed by m deviation vectors (2 < m < 2N) about any reference orbit and 
have been shown to: (a) Distinguish the regular or chaotic nature of the orbit faster than other methods, 
(b) identify the dimensionality of the space of regular motion and (c) predict the slow (chaotic) diffusion of 
orbits, long before it is observed in the actual oscillations. 

In the present paper, we improve the GALI method by introducing the Linear Dependence Indices 
(LDI m ). The new indices retain the advantages of thc GALI m and display the same values as GALI, in 
regular as well as chaotic cases. More importantly, however, the computation of thc LDI m is much faster 
in CPU time, especially if the dimensionality of phase space becomes large (N 3> 10). The main purpose 
of this paper, therefore, is to strongly advocate the use of LDI, for the most rapid and efficient study of the 
dynamics of multi - dimensional conservative systems. 

For the computation of thc LDI m we use information from the evolution of m > 2 deviation vectors 
from the reference orbit, as GALI does. However, while GALI requires the computation of many m x m 
determinants at every time step [27I 0] in order to evaluate the norm of the corresponding wedge product, 
LDI achieves the same purpose simply by applying Singular Value Decomposition (SVD) to the 2N x m 
matrix formed by the deviation vectors. LDI is then computed as the product of the corresponding singular 
values of the above matrix. This not only provides the same numerical values as the corresponding GALI m , 
it also requires much less CPU time. 

The paper is organized as follows: In section [2] we introduce the new index, explain in detail its compu- 
tation and justify its validity theoretically. In section [31 we demonstrate the usefulness of the LDI method, 
by applying it to the famous Fermi - Pasta - Ulam (FPU) lattice model of N dof, for small and large N. 
Finally, in section 0] we present our conclusions, highlighting especially the advantages of the new index. 



2 Definition of the Linear Dependence Index (LDI) 

Let us consider the 2N - dimensional phase space of a Hamiltonian system: 

H = H(q 1 (t),...,q N (t),p 1 (t),...,p N (t))=E (2) 

where qi(t), i = 1, . . . ,N are the canonical coordinates, Pi(t), i = 1, . . . , N are the corresponding conjugate 
momenta and E is the total energy. The time evolution of an orbit x{t) of §2§ associated with the initial 
condition: 

= (31 (*o), • • - ,QN(to),Pi(to), ■ ■ ■ ,2>jv(*o)) 
at initial time to is defined as the solution of the system of 2N first order differential equations (ODE): 

dqi{t) _ dH d Pi (t) _ dH 



dt dpi{ty dt dq l {t) 



' 1 V. (3) 



Eqs. ((3]) are known as Hamilton's equations of motion and the reference orbit under study is the solution 
x(t) which passes by the initial condition x(to). 

In order to define the Linear Dependence Index (LDI) we need to introduce the variational equations. 
These are the corresponding linearized equations of the ODE (J3J), about the reference orbit x(t) defined by 
the relation: 

dvi(t) 



dt 



J(x(t))-Vi(t), i=l,...,2N (4) 



2 



where J{x{t)) is the Jacobian of the right hand side of the system of ODEs calculated about the orbit 
x(t). Vectors Vi(t) = (ui,i(t), . . . , i>i 2Jv(*))i « = 1, • • • , 2N are known as deviation vectors and belong to the 
tangent space of the reference orbit at every time t. 

We then choose m € [2, 2N] initially linearly independent deviation vectors Vm(0) and integrate equation 
(U)) together with the equations of motion (|3J). These vectors form the columns of a 2N x m matrix A(t) 
and are taken to lie along the orthogonal axes of a unit ball in the tangent space of the orbit x(t) so that 
v m (0) are orthonormal. Thus, at every time step, we check the linear dependence of the deviation vectors 
by performing Singular Value Decomposition on A(t) decomposing it as follows: 

A(t) = U(t)-W(t)-V(t) T , (5) 

where U(t) is a 2N x m matrix, V(t) is an m x m matrix whose columns are the v m (t) deviation vectors 
and W(t) is a diagonal mxm matrix, whose entries wi(t), . . . ,w m (t) are zero or positive real numbers. 
They are called the singular values of A(t). Matrices U(t) and V(t) are orthogonal so that U T (t) ■ U(t) = 
V T (t) ■ V(t) = I, where I is the rectangular 2N x 2N unit matrix. 

We next define the generalized Linear Dependence Index of order m or LDI TO as the function: 

in 

LDI ro (i) = JI «;,•(*) (6) 
j=i 

with m = 2, 3, . . . , 2N, where N is the number of dof of 

The reason for defining LDI through relation © is the following: According to 2jj it is possible to 
determine whether an orbit is chaotic or lies on a d - dimensional torus by choosing m deviation vectors and 
computing the GALI,„ index. If GALI m ~ const, for m = 2, 3, . . . , d and for m > d decay by a power law, 
the motion lies on a d — dimensional torus. If, on the other hand, all GALI m indices decay exponentially 
the motion is chaotic. Thus, to characterize orbits we often have to compute GALI m indices for m as high 
as N or higher. 

A serious limitation appears, of course, in the case of Hamiltonian systems of large N, where GALIzv(t) 

involves the computation of = jwrfi determinants at every time step. For example, in a Hamiltonian 

system of N = 15 dof, GALI 15 (£) requires, for a given orbit, the computation of 155117520 determinants 
at every time step while LDI(t) = LDI 15 (i) requires only the application of the SVD method for a 30 x 15 
matrix A(t)l 

Clearly, at every point of the orbit x(t) the 2 < m < 2N deviation vectors span a subspace of the 2N - 
dimensional tangent space of the orbit, which is isomorphic to the Euclidean 2N - dimensional phase space 
of the Hamiltonian system |3]). Thus, if k of the m singular values tvk{t),k = 1, . . . , m are equal to zero, 
then k columns of matrix A(t) of deviation vectors are linearly dependent with the remaining ones and the 
subspace spanned by the column vectors of matrix A(t) is d(= m — k) - dimensional. 

From a more geometrical point of view, let us note that the m variational equations ^ combined with 
the equations of motion ([3]) describe the evolution of an initial m - dimensional unit ball into an m (or less) 
dimensional ellipsoid in the tangent space of the Hamiltonian flow. Now, the deviation vectors Vi{t) forming 
the columns of A(i) do not necessarily coincide with the ellipsoid's principal axes. On the other hand, in 
the case of a chaotic orbit, every generically chosen initial deviation vector has a component in the direction 
of the maximum (positive) Lyapunov exponent, so that all initial tangent vectors in the long run, will be 
aligned with the longest principal axis of the ellipsoid. The key idea behind the LDI method is to take 
advantage of this fact to overcome the costly calculation of the many determinants arising in the GALI m 
method and characterize a reference orbit as chaotic or not, via the trends of the stretching and shrinking 
of the m principal axes of the ellipsoid. 

Thus, LDI solves the problem of orbit characterization by finding new orthogonal axes for the ellipsoid 
at every time step and taking advantage of the SVD method. Since the matrix V in (|2|) is orthogonal, we 
have V T — V~ x , so that equation (|5|) gives: 

»4.2JVxm ' Vmxm = U^Nxm ' W mxm (7) 
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at every time step. Geometrically, Eq. {7J implies that the image formed by the column vectors of matrix V 

tVi ..... . 

is equal to an ellipsoid whose i principal axis direction in the tangent space of the reference orbit is given 

by: 

Wi ■ Ui (8) 

where Wi are the singular values of matrix A(t) and Uj is the i column of matrix U(t). This is, in fact, the 
content of a famous theorem stating that: 

Theorem 1 (f^J) Let A be a 2N x m matrix, and let U and W be matrices resulting from the SVD of A. 
Then, the columns of A span an ellipsoid whose i principal axis is u>; -Ui, where W = diag(wi, W2, ■ ■ ■ , w m ) 
(singular values) and {ui}™ =1 are the columns ofU. 

According to this theorem, the principal axes of the ellipsoid created by the time evolution of equation 
((4]) in the tangent space of the reference orbit x(t) at every time t, are stretched or shrunk, according to the 
singular values of Wi > 1 or Wi < 1 respectively for i = 1, . . . , m. 

If it so happens that k of the singular values Wi = as t grows, then the corresponding principal axes of 
the ellipsoid vanish and the ellipsoid is less than m - dimensional in the tangent space of the reference orbit 
because the corresponding deviation vectors of matrix A have become linearly dependent. 

Thus, two distinct cases exist depending on whether the reference orbit x(t) is chaotic or ordered: 

1. If the orbit is chaotic, the m deviation vectors become linearly dependent so that GALI m (i) — > 
exponentially [27| . Consequently, at least one of the singular values Wi(t),i = 2, . . . , m becomes zero 
and LDI m (i) = J]™ w 3 {t) -> (also LDI(t) -> 0) for all m > i. 



If the orbit is ordered (i.e. quasiperiodic) lying on a d - dimensional torus, there is no reason 23k l25| 
for the m deviation vectors to become linearly dependent, as long as m < d. No principal axis of the 
ellipsoid is eliminated, since all singular values Wi, i = l,...,m are nonzero and LDI m (t) fluctuates 
around nonzero positive values. On the other hand, for m > d, the singular values id,, i = d+ 1, . . . ,m 
tend to zero following a power law (27j . since m—d deviations will eventually become linearly dependent 
with those spanning the d - dimensional tangent space of the torus [13; Hi]. 



In the remainder of the paper, we apply the LDI indices and numerically demonstrate that: 

LDI m = GALI m , to = 2, . . . , 2N (9) 

for the same choice of m initially linearly independent deviation vectors fi(0), i = 1, . . . ,m. In particular, 
we present evidence that supports the validity of relation ^ and exploit it to identify rapidly and reliably 
ordered and chaotic orbits in a 1 - dimensional, N degree of freedom Fermi - Pasta - Ulam lattice under 
fixed and periodic boundary conditions 0; HJ. We propose that the validity of ([9]) is due to the fact that 
both quantities measure the volume of the same ellipsoid, the difference being that, in the case of the LDI, 
the principal axes of the ellipsoid are orthogonal. As we have not proved it, however, this is a point to which 
we intend to return in a future publication. 



3 Application to the FPU Hamiltonian System 

In this section, we apply the LDI method to the case of a multidimensional Hamiltonian system. Our aim is 
the comparison of its performance and effectiveness in distinguishing between ordered and chaotic behavior 
compared with Lyapunov exponents as well as the SALI and GALI methods. 

We shall use the TV dof Hamiltonian system of the ID lattice of the Fermi - Pasta - Ulam (FPU) [3 - 
model. The system is described by a Hamiltonian function containing quadratic and quartic nearest neighbor 
interactions: 

1 N N f\ 1 \ 

H * = 2 E ±2 3 + E ( 2 - + S x i +i - x i ]i = E (10) 

3 = l j=0 ^ > 
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Figure 1: The case of an ordered orbit: (a) The time evolution of the three maximal positive Lyapunov 
exponents, (b) The time evolution of GALI2 and LDI2. (c) The time evolution of GALI3, LDI3 and GALI4, 
LDI4. (d) The time evolution of GALI m . LDI m with m = 5, . . . , 8 and the corresponding slopes of the fall 
to zero. In all panels we have chosen a neighboring orbit at a distance ~ 2.1 from the OPM of the FPU 
Hamiltonian (|10|) with periodic boundary conditions for TV = 4 and E = 2. All axes are logarithmic. 



where Xj is the displacement of the j particle from its equilibrium position, ij is the corresponding 
conjugate momentum. /3 is a positive real constant and E is the constant energy of the system. 

We start by focusing on an ordered case choosing a neighboring orbit of the stable out of phase mode 
(opm) sum, which is a simple periodic orbit of the FPU Hamiltonian (fT0|) . This solution exists for 
every N, for fixed as well as periodic boundary condition (PBC): 

x N+1 (t)=x 1 (t),Vt (11) 

and is given by: 

Xj(t) = -x j+1 (t), Xj{t) = 0, j = 1, . . .,N,Vt. (12) 

In 0; Q the stability properties of the OPM mode with periodic boundary conditions were determined 
using Floquet theory and monodromy matrix analysis and the energy range < E(N) < i?^? PM (A) over 
which it is linearly stable was studied in detail. 
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It is known that for N = 4 and (3 = 1, the solution (TT2"]) with periodic boundary condition (TT2")) is 
destabilized for the first time at the critical energy 

£OPM K 4 51 _ Below thig critical 

energy, the OPM is 

linearly stable and is surrounded by a sizable island of stability. By contrast, for E > E® PM , the OPM is 
linearly unstable with no island of stability around it. 

In Fig. [IJa), we have calculated the three maximal Lyapunov exponents of a neighboring orbit located at 
distance ps 2.1 away from the OPM at E = 2 < E° PM . At this energy, the OPM is linearly stable and thus 
all Lyapunov Exponents tend to zero following a simple power law. Next, in Fig. [ljb), we compute GALI2 
and LDI2 for a final integration time t = 8 x 10 4 and observe that GALI2 and LDI2 practically coincide 
fluctuating around non zero values indicating the ordered nature of the orbit. GALI2 needs 558 seconds of 
computation time while LDI2 takes about 912 seconds in a Pentium 4 3.2GHz computer. 

In Fig. [IJc), we compute GALI3, LDI3 and GALI4, LDI4 for the same energy and initial condition. We 
see once more that GALI3, LDI3 and GALI4, LDI4 coincide fluctuating around non zero values. The GALI3 
computation now takes about 1044 seconds, LDI3 about 838 seconds, GALI4 needs 898 seconds and LDI4 
753 seconds. 

Finally, in Fig. |TJd) , we present GALI m , LDI m with m = 5, ... ,8 as a function of time for the same 
energy and initial condition. We observe again that GALI r „ and LDI m with m — 5, . . . , 8, have the same 
values and tend to zero following a power law of the form t~ 2 ( k ~ N ) , All these results are in accordance 
with the formulae reported in [27J and suggest that the torus on which the orbit lies is 4 - dimensional, as 
expected from the fact that the number of dof of the system is N = 4. 

In Q we also studied the stability properties of a different simple periodic orbit of FPU called the SPOl 
mode with fixed boundary conditions (FBC). Using monodromy matrix analysis we found that for TV = 5 and 
(3 = 1.04, the SPOl mode with FBC is destabilized for the first time at the critical energy E^ PQ1 ps 6.4932. 

Thus, in order to study a chaotic case where things are different, we choose initial condition at distance 
of ps 1.27 x 10 -4 from the SPOl orbit at the energy E = 11, where it is unstable. 

In Fig. EJa), we calculate Lyapunov exponents of the above mentioned orbit and find that the four 
maximal Lyapunov exponents tend to positive values. This is strong evidence that the nature of the orbit is 
chaotic. Next, in Fig. [^b) we calculate GALI2 and LDI2 up to t = 1200. We see the indices again coincide 
and tend to zero as oc e-^i-"^)* (solid straight line), as predicted by our theory 2(| 27 1. In this figure, we 
find tj\ ~ 0.124 and 02 ps 0.056 for time t = 71. The corresponding CPU time required for the calculation 
of all indices does not differ significantly, as they become quite small in magnitude, rather quickly. 

Nevertheless, LDI 2 requires less CPU time than GALI 2 . In Fig. [UJc), we calculate GALI 5 and LDI 5 
for the same energy and initial condition as in the previous panels. We observe now that GALI5 and 
LDI5 coincide falling to zero as GALI 5 cx e-K^i-^H^i-^H^i-^H^i-a.,)]* ^ solic j s t ra ight line) where 
u\ ps 0.197, cr 2 ~ 0.095, ct 3 ps 0.047, a 4 ps 0.026 and cr 5 ps 0.022 for t ps 44. Clearly, GALI 5 and LDI 5 
distinguish the chaotic character of the orbit faster than GALI2 or LDI2. This is so, because GALI2 or 
LDI 2 reaches the threshold 10" 8 [II [H; [27| for t « 750 while GALI 5 and LDI 5 for t « 35! The CPU times 



required for the calculation of GALI5 and LDI5 up to t — 80 are approximately 1.5 seconds each. 

Thus, we conclude from these results that the LDI method performs at least as well as the GALI, 
predicting correctly the ordered or chaotic nature of orbits in Hamiltonian systems for low dimensions, i.e. 
at 2, 4 and 5 degrees of freedom. However, in higher dimensional cases, GALI indices become very impractical 
as they demand the computations of millions of determinants at every time step making the LDI method 
much more useful. 

In order to show the advantages of the LDI method concerning the CPU time, we repeat the above 
analysis for the same Hamiltonian system (I10p . but now for N = 15 and energy E = 2, and for an initial 
condition very close to the unstable SPOl [4|. 

In 0] it has also been shown that for N = 15 and [3 = 1.04, the SPOl with fixed boundary conditions 
destabilizes at the critical energy E c ps 1.55. Thus, for energies smaller than E c , SPOl is linearly stable, 
while for E > E c it is unstable and is surrounded by a chaotic region. 

In Fig. [3Ja) we depict the time evolution of the five maximal Lyapunov exponents which converge to 
positive values for high enough t suggesting that the neighboring orbit is chaotic. In the second panel of the 
same figure we present the evolution of GALL3 and LDIs together with the approximate exponential law. 
We remark once more that the values of the corresponding indices coincide until they become numerically 
zero. More interestingly, the CPU time required for the calculation of GALIs up to t ps 100 is about 186 
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Figure 2: The case of a chaotic orbit: (a) The time evolution of the four maximal Lyapunov exponents, (b) 
The time evolution of GALI2, LDI2 follows the approximate formula e~( CTl ~ cr2 ) t where o\ « 0.124 (solid 
straight line) and 02 ~ 0.056 for £ = 71. (c) The time evolution of GALI5, LDI5 follows the approximate 
formula oc e -[(^-^)+^i-^)+^i-^)+{<n-^)]t _ e -o.069t ( golid straight line ) where CTl K .197, cr 2 « 0.095, 

03 ~ 0.047, 04 « 0.026 and 05 ps 0.022 for t w 44. We have used, in all figures, the same orbit of a distance 
of 1.27 x 10~ 4 from the SPOl of the FPU Hamiltonian system (|10[) with fixed boundary conditions for N = 5 
and E = 11. 

seconds while for the LDL3 it takes only one second! This difference is very important, showing why LDI is 
preferable compared to the corresponding GALI index in Hamiltonian systems of many degrees of freedom. 

4 Conclusions 

In this paper we have introduced a new method for distinguishing quickly and reliably between ordered 
and chaotic orbits of multidimensional Hamiltonian systems and argued about its validity justifying it in 
the ordered and chaotic case. It is based on the recently introduced theory of the Generalized Alignment 
Indices (GALI). Following this theory, the key point in the distinction between order and chaos is the linear 
dependence (or independence) of deviation vectors from a reference orbit. Consequently, the method of LDI 
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Figure 3: (a) The time evolution of the five positive Lyapunov exponents, (b) The time evolution of the 
GALI 8 , LDI 8 follows the approximate formula oc e -l( a i-^)+^i-^3)+-+(<^i-^a)}t ^ e -o.385t ( solic j straight 
line) where a Y w 0.061, a 2 « 0.011, <j 3 ~ 0.006, a 4 ~ 0.005, a 5 w 0.005, er 6 w 0.004, a 7 « 0.004 and 
as ~ 0.004 for time t ss 141. In all panels we have used initial conditions at a distance of 9 x 10 -5 from the 
SPOl orbit of Hamiltonian system (fT0|) with fixed boundary conditions, N — 15 and E = 2. 



takes advantage of this property and analyzes m deviation vectors using Singular Value Decomposition to 
decide whether the reference orbit is chaotic or ordered. If the orbit under consideration is chaotic then 
the deviation vectors are aligned with the direction of the maximal Lyapunov exponent and thus become 
linearly dependent. On the other hand, if the reference orbit is ordered then there is no unstable direction 
and to = 1,2, d < ./V deviation vectors are linearly independent. As a consequence, the LDI of order to 
(LDI m ) becomes cither zero if the reference orbit is chaotic or it fluctuates around non zero values if the 
orbit is ordered if to < d. 

After introducing the new method, we presented strong numerical evidence about its validity and effi- 
ciency in the interesting case of multidimensional Hamiltonian systems. One first main result is that GALI m 
and LDI TO coincide numerically for the same to number of deviation vectors and for the same reference orbit. 
Moreover, it follows that it is preferable to use the LDI method rather than the equivalent GALI method 
especially in the multidimensional case of Hamiltonian systems, since the LDI needs considerably less CPU 
time than the corresponding GALI method for the same number of deviation vectors. 
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